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We discuss the thermal conduction and convection of thermally relativistic fluids 
between two parallel walls under the gravitational force, both theoretically and nu¬ 
merically. In the theoretical discussion, we assume that the Lorentz contraction is 
ignored and spacetime is flat. For understanding of the thermal conduction and con¬ 
vection of thermally relativistic fluids between two parallel walls under the gravita¬ 
tional force, we solve the relativistic Boltzmann equation using the direct simulation 
Monte Carlo method. Numerical results indicate that strongly nonequilibrium states 
are formed in vicinities of two walls, which do not allow us to discuss the transition of 
the thermal conduction to the thermal convection of thermally relativistic fluids un¬ 
der the gravitational force in the framework of the relativistic Navier-Stokes-Fourier 
equation, when the flow-field is under the transition regime between the rarefied and 
continuum regimes, whereas such strongly nonequilibrium states are not formed in 
vicinities of two walls under the nonrelativistic limit. 
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I. INTRODUCTION 


In the early epoch of universe, the space is filled with high energy partons, whose ther¬ 
mal energy is larger than A QCD Q, which enables us to asstune the asyndetic fteedotn 
of partons. In this paper, we assume that partons are approximately expressed with hard 
spheres under the asymptotic freedom of partons in the early epoch of universe to simplify 
our discussion. Of course, the collisional differential cross section depends on the momentum 
transferred between two colliding partons, and the collisional deflection-angle also depends 
on the momentum, which is transferred between two colliding partons in the energy-regime 


of the asymptotic freedom 


1] [ 12 ]. In particular, we characterize the thermal energy of par- 


tons using the thermally relativistic measure, namely, y = me 2 / ( kd ) (m: mass, c: speed of 
light, 9: temperature, k: Boltzmann constant), whereas we assume that m is large enough 
to realize y —>■ 0 by not m —> 0 but 9 —> oo. Additionally, A QCD <C me 2 is assumed. Conse¬ 
quently, thermally relativistic fluids are characterized by y < 100. Additionally, we assume 
that the asymptotic freedom of partons always consists under y < 100. 

In our previous studies, we discussed thermally relativistic fluids with dissipation, which are 
composed of hard spherical partons, such as the shock layer or thermal fluctuation under 
the equilibrium state [sj] . Meanwhile, the characteristics of the thermal conduction and con¬ 
vection of thermally relativistic fluids under the gravitational force, which are composed of 
hard spherical partons, have not been understood, adequately. In this paper, we investigate 
the thermal conduction and convection of thermally relativistic fluids under the gravita¬ 
tional force in flat spacetime for understanding of the thermal conduction and convection of 
thermally relativistic fluids under the gravitational force in the early epoch of the universe. 
Of course, such an assumption of flat spacetime is inadequate to discuss the thermal con¬ 
duction and convection of thermally relativistic fluids under the gravitational force in the 
early epoch of universe, because spacetime is markedly curved in the early epoch of universe 
4], We, however, consider that our present study of the thermal conduction and convection 
of thermally relativistic fluids under the gravitational force in flat spacetime will be useful 
for understanding of the thermal conduction and convection of thermally relativistic fluids 
under curved spacetime in the early epoch of universe. 

The thermal conduction and convection of thermally relativistic fluids under the gravi¬ 
tational force in flat spacetime are studied by investigating the thermal conduction and 
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convection between two parallel walls, which exist at x = 0 and L, as shown in Fig. 1. 
Temperatures of two walls (0 V temperature of the cold wall at x\ = 0, O'?: temperature of 
the hot wall at x\ = L) are different, namely, 9™ < 0™, whereas the gravitational force is 
defined by g = (gi, < 72 , <73 ) = (<7,0,0) (0 < g). In particular, the thermal convection via such 
a gravitational force is called as Rayleigh Benard convection (fj . For example, two parallel 
walls correspond to thermal bathes (energy sources) in the early epoch of universe. 

In this paper, we investigate the characteristics of the thermal conduction and thermal 
convection, namely, Rayleigh Benard convection of thermally relativistic fluids, both the¬ 
oretically and numerically, whereas effects of Lorentz contraction are ignorable owing to 
|it| <C c ( u: flow velocity). Our previous studies [si indicated that the relativistic Navier- 
Stokes-Fourier (NSF) equation is efficient for the analytical comprehension of the weakly 
nonequilibrium state such as thermally fluctuations under the static state, when the Lorentz 
contraction is ignorable and thermally relativistic effects are dominant. 

This paper is organized as follows. In Sec. II, we define the Rayleigh number for thermally 
relativistic fluids on the basis of the relativistic NSF equation to introduce the theoreti¬ 
cal value of the critical Rayleigh number in appendix A. In Sec. II, temperatures on two 
walls, which define the Rayleigh number, are not defined, analytically. Thus, we consider 
temperature jumps on two walls on the basis of Grad’s 14 moments equation coupled with 
Fourier law in Sec. III. Here, we must notice that temperature jumps on two walls are cal¬ 
culated using the method 6], which has been used to calculate the temperature jump of the 
non-relativistic gas on the wall as a half-space problem [lj, whereas Israel-Stewart equation, 
which was applied to describe the boundless Truesdell’s uniform shear flow by Yano-Suzuki 
sj], cannot be applied to describe the profiles of Grad’s 14 moments on the wall, because 
the distribution function is discontinuous on the wall. In Sec. Ill-A, we find analytical 
solutions of the density and temperature under the pure conductive state under two limits, 
namely, nonrelativistic limit (y y °°) and thermally relativistic limit (y —>■ 0). Finally, we 
calculate temperature jumps on two walls in the framework of Grad’s 14 moments coupled 
with the Fourier law under two limits using the analytical solution of the temperature un¬ 
der the pure conductive state. I 11 Sec. IV, we numerically investigate the characteristics 
of the thermal conduction and convection of thermally relativistic fluids under the gravita¬ 
tional force by solving the relativistic Boltzmann equation on the basis of direct simulation 
Monte Carlo (DSMC) method 9] ( 3 ] . Numerical results, however, indicate that the tran- 
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sition of the thermal conduction to the thermal convection of thermally relativistic fluids 
under the gravitational force cannot be discussed in the framework of the relativistic NSF 
equation owing to strongly nonequilibrium states in vicinities of cold and hot walls, when 
the flow-field is under the transition regime between the rarefied and continuum regimes. 
In particular, the temperature jumps on two walls increase in accordance with the decrease 
in the temperature, because the representative relaxation rate [loj] decreases in accordance 
with the decrease in the temperature, when Knudsen number is fixed. As a result, we cannot 
use the Rayleigh number as a parameter, which characterizes the transition of the thermal 
conduction to the thermal convection, under the transition regime between the rarefied and 
continuum regimes. Finally, we make concluding remarks in Sec. V. 


II. DEFINITION OF RAYLEIGH NUMBER FOR THERMALLY 

RELATIVISTIC FLUIDS 


The relativistic NSF equation is written for thermally relativistic fluids with zero Lorentz 
contraction in acausal form as 

dp 


dt 

dpuiG 
dt 
dd 
dt 


= - V • (pu ), 
d 


dx k 


(p$ik + pUiU k G + II ik ) + pg h 


pc v 


= — pc v u • V0 — V • q — pV • u, 


(1) 

( 2 ) 
( 3 ) 


where p is the density, u = (wi,« 2 , W 3 ) is the flow velocity, p is the static pressure, g l is the 
gravitational components, G = K 3 (y) /K- 2 (y), in which K n is the n-th order modified Bessel 
function of the second kind, II,is the deviatoric stress tensor, q = ( 91 , 92 , 93 ) is the heat 
flux vector, and c v = k (y 2 + 5Gy — G 2 y 2 — 1) is the specific heat at the constant volume. 
The derivation of Eqs. (l)-(3) is described in appendix C in detail. The formal difference 
between the thermally relativistic NSF equation and nonrelativistic NSF equation is a term 
G in Eq. (2). Consequently, the form of thermally relativistic NSF equation coincides with 
the form of the nonrelativistic NSF equation, when G — 1. 

Tlij in Eq. (2) and q in Eq. (3) are written in the framework of the NSF law as follows 


n ij = -rj 


dui duA ( 2 \ dui 

&‘ + 3^j"( C "3’T«5V 


q = -AV0, 


( 4 ) 

( 5 ) 
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where g is the viscosity coefficient, ( is the bulk viscosity, and A is the thermal conductivity, 


which are calculated for hard spherical partons by Cercignani and Krerner 
of Israel-Stewart theory. 

Provided that the flow is incompressible, we can rewrite Eqs. (l)-(3) as 

V • u — 0, 

duiG 0 2 

p—— = - vp - pu k —UiG + r/V Ui + pgi, 

09 


10] on the basis 


pc v — = — pc v u • Vd + AV 2 0 — p'V ■ u , 
Ot 


( 6 ) 

(7) 

( 8 ) 


Next, we assume that quantities can be expressed with those under the background state 
and their perturbations as 


p = p 0 + Sp, u = u 0 + Su, p = p 0 + Sp, 

9 = 9 0 + 89, G = G 0 + SG ~ G 0 - G 0 ^, 

9 o 

where 9 0 = 6L and p 0 = PqROq (R — k/m), in which 9_ = 9\ X1=0 . 

Equation (7) yields the following relation under the static state, namely, u = 0 

Pod = Vp 0 , 


(9) 


( 10 ) 


Subtracting Eq. (10) from Eq. (7) and dividing both sides of Eq. (7) by poGo, we rewrite 
Eq. (7) using Eq. (9) as 

OSui 


(1 + Sp) 


Ot 


0 


= -V5p - (1 + Sp) Suk^—Sui + rjX7 2 Sui + 8pep 

dxk 


( 11 ) 


where we used approximations Sui ~ 0 an d SpSGgt ~ 0, and relations Sp = 

Sp/po, SG = SG/Gq, Sp = Sp/ {poG 0 ), 77 = 77 / (p 0 Go), and g t = gi/G 0 . 

Setting Sp ~ 0 in the first term of the left hand side of Eq. (11) and second term of the 

right hand side of Eq. (11), Eq. (11) is rewritten as 


OSUi 0 _ r r~~ 

-yj- = -vSp - 8u k —8ui + ri\7~8ui + 5pg u 


( 12 ) 


Substituting Eq. (9) into Eq. (8), we obtain the linear balance equation of energy by 
neglecting nonlinear terms 

086 


Ot 


= A V 2 56 - 5u ■ V89, 


(13) 
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where we used the relation — p'V ■ v = —p {c p — c v ) + u ■ VP) and A = A/ (c p p 0 ). 
To nondimensionalize Eqs. (12) and (13), we introduce following relations 

5p = ^ 2 ^5p, 5u = $/5u, 56 = (9 1 — 9 0 ) 59, x = Lx, t = —t, 

where 9\ — 9 + , in which 6 + = 9\ X1=L . 

Substituting Eq. (14) into Eq. (12), we obtain using V = TV as 

W 2 d5Ui W 2 d ^ 2 Pl’- -,=72 r- r~~ 

~P -aT =- —5u k ——5ui - — — V5p + — r?V 5ui + 5pg h 

L ot L uXk L Go L/ 

where Pr = pc p /X is Prandtl number. 

Here, we use Boussinesq’s approximation: 


(14) 


(15) 


P — Po {l + a ($1 — Po) ) 


(16) 


where a = —l/p 0 (8p/89) p=po . 
From Eq. (16), we readily obtain 


5p = -—— = a (0i - P 0 ) 59. 
Po 


(17) 


Substituting Eq. (17) into Eq. (15) and dividing both sides of Eq. (15) by W 2 /L, we obtain 


- ^V5p + V 2 ^ + A a ( 01 _ p 0 ) 59~ 9i , 

ot oXk Go °u L °// A 

Here, we rewrite Eq. (18) using the relation % = X/L as 

c)&u ■ c) /s _ /s _ /v — 

= —5uk ——— PrVhp' + PrV 2 (7uj + Pr —ck (Pi — P Q ) 59 gi , 


dt 


dxi 


Xp 


(18) 


(19) 


where 5p’ = 5p/po and Pr = Pr /Gq. 

Dividing both sides of Eq. (19) by Pr, we obtain 


1 85Ui 1 _ 8 - c- ^ c-77 

——-=- = — —5u k -^—5ui — W5p + \7~5ui + Ra APe*, 
Pr dt p r 8x k 


( 20 ) 

where Ra = g L 3 a (Pi — Po) / (^Ar)J is the Rayleigh number for thermally relativistic fluids, 
and e is the unit vector. 

Finally, Ra for thermally relativistic fluids is 

g L 3 p 2 c v a (Pi - P 0 ) 


Ra = 


77 A 


( 21 ) 


where c p = k (y 2 + 5 G\ — G 2 \ 2 ) is the specific heat at constant pressure. 

Provided that thermally relativistic fluids are composed of hard spherical partons with mass 
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m and diameter d, we obtain Ra using definitions of g and A for hard spherical partons, 


which were calculated by Cercignani and Krerner on the basis of Israel-Stewart theory 
as 


loj, 


Ra 


256 ^ , 

45 FrKn 2 V ' ’ 


(r = 0 O M) 


{(X 2 + 2) K 2 (2 X ) + 5 X K 3 (2 X )} {(15 X 2 + 2) K 2 (2 X ) + x (3y 2 + 49) K 3 {2 X )} 


X 7 K 3 ( X y ( X K 2 ( X y + 5K 3 ( X )K 2 ( X ) - X K 3 ( X ) 2 ) 
0.29 (x —>• 0) 

0.38 (x —> oo), 


( 22 ) 


where Fr = (2 kQi/m) / ( gL) is Froude number, Kn = 1/ (y/2no7Td 2 L) is Knudsen number, 
in which no = po/m is the number density, and a = 1/9. Stefanov et al. calculated the 
Rayleigh number for nonrelativistic fluids as Ra = 2.16 (1 — r) / (FrKn 2 ) j^J, which coincides 
with Ra under X —> oo in Eq. (22). 

In this paper, we restrict ourselves to the two dimensional flow. Thus, we consider the flow 
on (x\ , ,x' 2 )-plane. Additionally, the external force (g,) acts along the positive aq-direction, 
namely, g 2 = g 3 = 0 A 0 < g 3 . Forms of balance equations of the momentum and energy 
in Eqs. (20) and (13) are equivalent to those for nonrelativistic fluids, when we set Pr = Pr 
in Eq. (20) and A = A in Eq. (13). Therefore, the calculation of the critical Rayleigh 
number, which indicates the transition of the thermal conduction to the thermal convection, 
is obtained in a similar way, which is used for the calculation of the critical Rayleigh number 
of nonrelativistic fluids, as shown in appendix A. Actually, Stefanov et al. jf3] indicated 
that the critical Rayleigh number, which is numerically obtained using the DSMC method, 
coincides with the theoretical value 1708 in Eq. (A26) when Kn = 0.01. Finally, the critical 
Rayleigh number of thermally relativistic fluids for three types of the boundary condition, 
namely, rigid (hot wall)-rigid (cold wall), slip (hot wall)-rigid (cold wall), and slip (hot wall)- 
slip (cold wall), are quite same as those of nonrelativistic fluids, as discussed in appendix A. 
In above discussion, 6± are not defined. Then, we consider temperature jumps on two walls 
by calculating analytical solutions of the density and temperature under the pure conductive 
state in the next section. 





A Analytical solutions of density and temperature under pure conductive state 


III. TEMPERATURE JUMPS ON TWO WALLS IN THE FRAMEWORK OF 
GRAD’S 14 MOMENTS EQUATION COUPLED WITH FOURIER LAW 

In this section, we calculate temperature jumps on two walls in the framework of Grad’s 
14 moments equation coupled with the Fourier law using analytical solutions of the density 
and temperature under the pure conductive state, where we discriminate the pure conduc¬ 
tive state from the conductive state (u ^ 0) by u = 0. 


A. Analytical solutions of density and temperature under pure conductive state 


Substituting u = 0 into Eqs. (2)-(3), we readily obtain following two relations using Eqs. 
(4) and (5) and A for hard spherical partons on the basis of Israel-Stewart theory 10] under 
the steady state. 

1 dp/x{x i) 2 


p(x i) dx\ 


Fr’ 


(23) 


d 
dx i 
= 0, 


X (u) {x(xi)K 2 (x (ah)) 2 + 5 K 3 (y(x x )) K 2 (y(xi)) - x (u) I<3 (x(xi)) 2 } 2 dx(x i) 1 


K 2 (x(xi)) 2 (5 X (xi)K 3 (2 X {xi)) + (x(xi) 2 + 2) K 2 (2 X (x 1 ))) 


dx i 


(24) 


where p = p/poo and 9 = 9/9^ are uniform in x 2 direction, in which p^ and 9^ are 
representative values of the density and temperature. As a result, p and 9 are functions of 
x\. 

Calculations of solutions of p(x i) and y (xi) i 11 Eqs. (23) and (24) are markedly difficult. 
Thus, we calculate solutions of p (xi) and x (^i) under the thermally relativistic limit (y —> 0) 
and nonrelativistic limit (y —y oo) in Eqs. (23) and (24). 

From Eqs. (23) and (24), we obtain solutions of p(x i) and 9 (x i) under the thermally 
relativistic limit (y —y 0) as 


P (u) = p{x i) = [-2c (axi + b)}~ 

9 (xi) = ax i + 6, 


(25) 

(26) 


where a = 9(1)—9(0), b = #(0), c = —Fr/2, and = J 0 ' p (xi) dx\ = (—2) «c (be )« — (ca + cb) i /2 
9(0)(= 9_) and 0(1)(= 9 + ) must be determined by considering the temperature jump on 









B Temperature jumps on two walls in the framework of Grad’s 14 moments equation coupled with Fourier 1 


the wall. 

From Eqs. (23) and (24), we obtain solutions of p(x i) and 6 (aq) under the nonrelativistic 

limit (x ~t o°) as 

P(%i) = , P(xi)=exp 

6 (xi) = (a'xi + b'Y , (28) 

where a' = 0( 1)§ — 0(0)§, b' = 0( 0)§, d = — 2/Fr, and c &' = f* p (xi) dx i. 


2 3r' i 

- In |a'x! + b'\ + —— (a'xi + b')* 

3 a' 


(27) 


B. Temperature jumps on two walls in the framework of Grad’s 14 moments 

equation coupled with Fourier law 


The temperature jump on the wall is explained in the framework of Grad’s 13 moments 
under the nonrelativistic limit with the good accuracy, when Kn < 0.01 jg]. Here, we 
consider temperature jumps on two walls in the framework of Grad’s 14 moments equation 
coupled with Fourier law using the analytical solutions of the density and temperature in 
Eqs. (25)-(28). 

At first, the distribution function of partons is approximated by Grad’s 14 moments equation 
(/ 14 ) on the basis of Israel-Stewart theory, which is defined by flo| 


/l4 = / (0) h + 


n 


1-5 g x -x 2 + g\ 2 


X 


p 20 G + 3y - 13 G 2 X - 2 G X 2 + 2 G\ 2 
15 G + 2x - 6 G 2 x + 5Gx 2 + X 3 - G 2 x 3 


1-5 g x -x 2 + g*x 2 
3y 6 G + x- G 2 X JT a X 

me 2 1 — 5 Gx — X 2 + G 2 x 2 m 2 c 4 


Qo 


X 


p X + 5G- G 2 x 


G 


1 


me* 


-,P 


m 2 c A 


Upp 0 !/ 6 


U a Upp a p 0 
77 (*P) X 


+ 


1 


p 2 G m 2 c 2 


p p 


(29) 


where p a is the four momentum of a parton, U a is the four velocity, n is the dynamic 
pressure, upp is the pressure deviator, and is the Maxwell-Jiittner function, which is 

dehned by 

n 


/ (0) = 


-u a p a 

-e ke 


(30) 


4'Km 2 ckdK 2 (y) 
where n = p/m is the number density. 

Equation (29) is reduced to the following form using relations, \u\ /c <C 1 and 7 r^)/ (pc 2 ) 1 
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(i 7 - j), which is numerically confirmed from our calculation, 

n+ x± - 5G ±x± - x± + G ±x± 


At = / (0) 


i + 


x 


n± 20 G± + 3x± — 13Gj_y± — 2 G±x± + 2G±X± 
15 G± + 2x± — 6 G\_x± + 5 G±x± + X± ~ G ±X± 

1 — 5 Gx± — X± + G ±X± 

1 &X±G± + 3x± - 3 G 2 ± x± 


i-5 G ±x ±-xi + Gixi 


■p + x{f> 


0\2 


<h 


x± 


n± X± + 5 G± - G 2 ±x± 


7T^~ -\/2 

(n ~l ~1 ~o\ , (**) a.± / ~i\2 

( G±p -p P ) + __ (P ) 


(31) 


where x± = me 2 / (fc0±), II = 11/ (n^mc 2 ), 7f {Q(3) = 7r (a/3) / (n^mc 2 ), p = p/ (n^mc 2 ), q a = 
q a / (noomc 3 ), h = n/rioo, and = p Q /(mc), in which the quantity with a subscript or 
superscript ± indicates the quantity on the cold wall (—) or hot wall (+). In Eq. (31), we 
used relations under the pure conductive state, namely, U a = (c, 0, 0, 0) and q a = (0, q 1 , 0, 0). 
Distribution functions of partons, which are reflected on perfectly diffusive hot and cold walls, 
(/±) are written as 


rw 
■ — 


n± 


4irm 2 ck8±K2 (x 


_ e -pox± ( p i < o) , 


(32) 


where x± = me 2 / (/c 0 ±) (9™: temperature of the cold wall, 9™: temperature of the hot wall). 
Continuum conditions of heat fluxes on two walls yield relations jg] 


i / i o/± d 3 p 
q± = c I pp f u ~r = c 


,1„0 *± d 3 P 
D(<,>)p 1 P 


P l p {) fu— + c 


V(<>>) 0 


1 0 r W d 3 p 

P P /±-n-, 

P U 


where we used T 01 ~ q 1 under \u \ /c 1 (T" /3 : energy-momentum tensor) 10 ]. 
Substituting /l in Eq. (31) and f± in Eq. (32) into Eq. (33), we obtain 


?± = ± 

+ 2 e _x± 


e x± (x± + 3x± + 3) 7T(ii)e x± (x± + &X± + 15x± + 15) 


(33) 


X±K 2 (x±) 2 x±K 3 (x±) 

n { — (x± + 3x± + &x± + 6x±)G± 

+ (x± — 4) (x± + f>X± + 15x± + 15) + (2x± + 15x1 + 39x± + 39y±) G±} 
{^2 (x±) [(4x| - 40 X |) G± - 4x 5 G| - 6x1 + 26x1^1] }" 


-l 


„ e _x ± {(x±) 2 + 3x± + 3} 

77 ---— 

(x±) 3 K-2 (x±) 


( 34 ) 
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Substituting the NSF law in Eqs. (4) and (5), and u\/c ~ 0 into Eq. (34), we obtain 


= ± 


dx i 


e-^( X | + 3 X ± + 3) _ e-^((xl) 2 +3x^ + 3) 

n± - \ , --- n± - 


X±K 2 (x= 


(x±) 3 K 2 (xl 


^ ( X± + 5 G± - G 2 ± x±) 2 X±K 2 (; x±f dx ± 1 

16 1Vii (x± + 2) K 2 (2x±) + 5x±iF 3 (2x±) dx i 


«-)• —— Kn 


± 


_ e" x± (xi + 3x± + 3) _ e 

n ± - o — 7 r-n± — 


-x± 


+ 3x^ + 3 


X±^2 (X: 


:± 


(x±) 3 -^2 (x±) 


(35) 


where A ± = X ± / (rioock). 

x(^i) and n(x i) can be calculated under two limits, namely, x ► 0 and x -*> oo in Eqs. 
(25)-(28). 

Substituting dx^/dxi = X+ 1 — X- 1 and n ± — Frxi 1 /^, which are obtained under the 
thermally relativistic limit (x —>■ 0) from Eq. (26), into Eq. (35), we obtain 

__3_ k (x± + 5 G± - G\x±) 2 X±K 2 (x±) 2 / _i _ _n 
16 11 (xi + 2)/l 2 (2x ± ) + 5 x±/l 3 (2x±) U+ ’ 


± 


;- M »- as * 11 - vw'"- " - ' 




X±K 2 (x 


± 


(xl) K 2 (xt) 


(36) 


where Fr = 2/ (gx+) (g = c 2 /L). 

Equation (36) defines x± f° r the pure conductive state under the thermally relativistic limit 
(x —$■ 0). Similarly, we obtain the equation for the pure conductive state from Eqs. (27), 
(28) and (35), which defines x± under the nonrelativistic limit (x —>■ oo) as 


1 Vn (x± + 5 G± - G 2 ± x±) 2 X±K'2 ( X± f v 1 
8 n (x± + 2) K 2 (2x±) + 5 x±K 3 (2x±) ° f/a'&± + V 


± 


. e_x± (x± + 3x± + 3) _ e~ x ± {(x±) 2 + 3 X ± + 3} 

n± x 3 ±k 2 (x±) n± (xZ) 3 k 2 (xX) 


(37) 


where h± = p±/m is obtained by Eq. (27), respectively, and © + = 1 and ©_ = 0. 

We investigate two parameters, which evaluate temperature jumps on two walls, namely, 
A 9 + = 0+/9™ — 1 = X+/X+ ~ 1 and A6L = 1 — 6L/6 1 ™ = 1 — X-/X- on two walls by setting 
X+/X- = 0-1 {6+/0™ = 10) and v^Kai = 0.01. The left frame of Fig. 2 shows A 6 1 } (yi axis) 
versus xZ (ah axis), and A 9]^ (y 2 axis) versus x+ (x 2 axis) in the range of 10~ 3 < X- A 10 _1 , 
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which are obtained by Eq. (36), when g = 10“ 4 , 100 and 1000, where the superscript 14 
means that the value is obtained using Grad’s 14 moments coupled with Fourier law. As 
shown in the left frame of Fig. 2, A# 1 : 1 decreases, as x~ (&-) decreases (increases), when 
g = 10~ 4 , 100 and 1000, whereas dAd^/dx- decreases, g increases. A#/ 4 slightly increases, 
X- ($-) decreases (increases), when g = 10 -4 and 100, whereas A#;) 4 slightly decreases, x~ 
(0™) increases (decreases), when g = 1000. The right frame of Fig. 2 shows — Ad 14 (y 1 axis) 
versus X- ( x i axis), and — Ad+ (y 2 axis) versus x+ ( x 2 axis) in the range of 20 < X- A 150, 
which are obtained by Eq. (37), when g = 10” 4 , 10~ 3 and 0.01. As shown in the right frame 
of Fig. 2, — Ad] 4 decreases, as x~ (&-) increases (decreases), when g = 10~ 4 , 10~ 3 and 
10~ 2 , whereas —Ad^/dx- decreases, as g decreases. Additionally, —A0/_ 4 decreases, as x+ 
(d™) increases (decreases), when g = 10 -4 and 10~ 3 . —A d 1 ^ decreases, as x+ ($+) increases 
(decreases) in the range of 2 < x+ A 6 , when g = 10 ~ 2 , whereas —A d]^ increases, as x+ 
(d™) increases (decreases) in the range of 6 < x+ A 15, when g = 10~ 2 . 


IV. NUMERICAL STUDY OF THERMAL CONDUCTION AND 
CONVECTION OF THERMALLY RELATIVISTIC FLUIDS BETWEEN TWO 
PARALLEL WALLS UNDER GRAVITATIONAL FORCE 


In this section, we numerically analyze the thermal conduction and convection of ther¬ 
mally relativistic fluids under the gravitational force between two parallel walls by solving 
the relativistic Boltzmann equation on the basis of the DSMC method. The details of the 
DSMC method are described in appendix B. As the numerical condition, we consider the 
rectangular domain, which is stretched by (xi,x 2 ) = ^L,4L^, as shown in Fig. 1, in which 
L = L/Loo = 1. 128 x 128 grids are equally spaced in the rectangular domain. As a result, 
the cell size is characterized by Aaq = 1/128 and Ax 2 = 1/32. 1638400 sample particles 
are set in the rectangular domain. Additionally, -\/2Km =1/ (■ nd 2 n 00 L 00 ) = 0.01 indicates 
that flow-field is under the transition regime between the rarefied and continuum regimes, 
where d is the diameter of a parton. Meanwhile, we must notice that the representative 
relaxational time f = y^Km/ ( v s /c ) (v s : speed of sound) 10] increases in accordance with 
the increase of y, because v s decreases in accordance with the increase of y. Therefore, 
temperature jumps on two walls depend on both of x± an d Kn. The periodic boundary 
condition is applied to x 2 = 0 and x 2 = 4, whereas partons are reflected on two perfectly 
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diffusive walls at X\ — 0 and x\ = 1. All the numerical tests are shown in Tab. 1. As nu¬ 
merical tests, we consider six cases, namely, case (A) (x+,X-) = (0.00029,0.0029), case (B) 
feX-) = (0.0029,0.029), case (C) = (0.01,0.1), case (D) (y^y™) = (0.25,2.11), 

case (E) (x+,X-) — (2.11,21.1) and case (F) (x+, X—) — (10,100.29), as shown in Tab. 1. 
Each of cases includes four tests, (i)-(iv). Test (i) corresponds to the conductive state (C), 
test (ii) corresponds to the weakly wavy state (WW), test (iii) corresponds to the obscure 
formation of four vortices (4 — V 1 ), and test (iv) corresponds to the clear formation of four 
vortices (4 — V 2 ). In the next subsection, we discuss conductive states in tests (i) of cases 
(A)-(F). 


A. Numerical results of conductive state 

Table 1 shows A 6 L and A 9 + together with A 6 U and A A) 4 , which are calculated using 
Eq. (36) for cases (A)-(C) and Eq. (37) for cases (D)-(F), where we calculated A 6 L and 
A 9 + by setting X- — X (Aaq) and x+ — X (1 — Aaq), because we cannot calculate tempera¬ 
tures on two walls owing to the characteristics of the DSMC method. As shown in Tab. 1, 
A 6 L ~ 0.06 in tests (i) of cases (A)-(C), whereas Ad 14 increases, as x~ increases in tests (i) 
of cases (A)-(C), as discussed in Sec. III-B. A 6 + increases in tests (i) of cases (A)-(C), as 
y+ increases, whereas A9]^ slightly increases in tests (i) of cases (A)-(C), as x+ decreases. 
The orders of A 6 L and A 6 + are, however, similar to those of A 6 14 and A 6 l + in tests (i) of 
cases (A)-(C). Additionally, A 6_ and A 6 + are similar to A9 14 and A 6]^ in test (i) of case 
(D). Meanwhile, A 6 L <A A 8 14 and A6]^ <C A 9 + are obtained in tests (i) of cases (E) and 
(F). Snch marked differences between A 6_ and AfA 1 or A 6 + and A A ) 4 might be explained 
by two reasons. One is that the representative relaxational time f markedly increases in 
accordance with the increase of y owing to the decrease of v s , when 1 < y. As a result, the 
nonequilibrium effect beyond Grad’s 14 moment approximation markedly increases when 
y increases from y = 1, even when Kn corresponds to the transition regime between the 
rarefied and continuum regimes. The other reason for such marked differences between A 6_ 
and A 6 14 or A 6 + and A 6 l + is that continuum conditions of q 1 on two walls and zero-flow 
condition under the pure conductive state are insufficient. 

Figure 3 shows u\ (yi axis) and n (y 2 axis) versus 1 — x\ in tests (i) of cases (A)-(F), where 
u = u/c. As shown in Fig. 3, the flow toward the hot wall (1 — x\ = 0) is markedly 
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accelerated in the vicinity of the hot wall in tests (i) of all the cases (A)-(F). Similarly, the 
flow toward the cold wall (1 — x± = 1) is also markedly accelerated in the vicinity of the 
cold wall in tests (i) of all the cases (A)-(F). In particular, the flow toward the hot wall is 
a reverse of the heat flow toward the cold wall. Such a marked acceleration in vicinities of 
two wall presumably yield the strongly nonequilibrimn state in the vicinity wall, which does 
not allow us to assume that the temperature jump can be described by the continuum con¬ 
dition of q 1 on two walls under the pure conductive state on the basis of Grad’s 14 moments 
coupled with Fourier law. II decreases in vicinities of two walls in cases (D)-(F), whereas II 
decreases in the vicinity of the hot wall, and markedly increases in the vicinity of the cold 
wall. Finally, II ^ 1 is obtained in all the cases. 

Figure 4 shows 7Tpq and q 1 versus 1 — x\ together with 7fA S A and g 4 SF , which are calculated 
by substituting numerical results of U\ and 6 into Eqs. (4) and (5). As shown in Fig. 4, 
< 17T(ii) | is obtained in vicinities of two walls, whereas ffpi) ~ 0 is obtained in the 


7T 


(U), 

regime except for vicinities of two walls. 17r(n) | / in the vicinity of the hot wall in¬ 

creases, as x+ increases, |g 4 | < |g 4 SF | is obtained in the vicinity of the hot wall, whereas 
I^nsfI < I? 1 ) is obtained in the vicinity of the cold wall. From above comparisons of 7T(n) and 
q 1 with 7f^ and g 4 SF , we confirm that strongly nonequilibrium states beyond the relativistic 
NSF equation are formed in vicinities of two walls in tests (i) of cases (A)-(F). 

Here, we must confirm whether marked differences between A 9_ and Aid 4 or A 9 + and A0+ 4 
in cases (E) and (F) are dismissed by considering no permeation of the flow in normal direc¬ 
tions to two walls ((hi)^ = 0) and continuum conditions of 7r^ together with continuum 
conditions of q± on two walls. 

In cases (E) and (F), we can assume that the number of partons, which satisfy the relation 
Xf^Ui <C 1, is much larger than the number of partons, which do not satisfy the relation 
XP 1 ^i <!• As a result, we can obtain following form from Eq. (31) by setting n ~ 0 


fu = fs% (! + X±P 1 uf) 


1 + 


9 ? 


X± 


n± X± + 5 G± - G 2 ± x± 


{G±p 1 -p 1 (P° + ufp 1 )} 


rfn) X± 

n± 4G± 


(?y 


( 38 ) 
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where we used relations, U a = U a /c ~ (1, iq, 0, 0), II ~ 0 and 

f± = fs°i exp {xP'uf) ~ fs% (! + X±P l uf) , 


AO) _ 


n± 


c(0) 

S,± 

-°P0 

k6± 


A'Rm 2 ck6±K 2 (x=t) 

Zero-flow in normal directions to two walls yields relations 

d 3 p 


0 = c f P X ft —+ c f P { f±^~- 


(39) 


(40) 


Continuum conditions of 7T(n) on two walls yield relations 


ttJ 1 ) = c / PiPif: 


± d 3 p 

14 ^T 


= c 


' 0(<,>)p 1 


PiPif: 


± d 3 p 


+ c 


V(<,>)0 


rw d3 P 
PlPlf±— Q-, 

pO 


(41) 


where we used Tu ~ 7T(n) under \u\ /c 1. 
Substituting ft in Eq. (38) into Eq. (33), we obtain 

e ~ X± ^fn) {2e x± x±UiK 4 (x±) ± Ct) 


& = 


± 


2X^3 (x±) 

h±e~ x± Cf fi±e~ x "Cf u 


+ n±G±Ui 


9± 


B 


(x±) 3 k 2 (x±) {x w ±?K 2 {xl) 

I {3 (uf) 2 xilu(x±) ± e~ x± uf (4 Ct - ±X±CfG± + 2 Cf) } , 


(42) 


where Bf = x±-Mx±) (x± + 5 G± - x±&±), c i = X± + 6x± + 15x± + 15, Cf = x± + 
9x 3 ± + 39x± + 90x4 + 90 and C± = xl + 3x± + 3. 

Substituting ft in Eq. (38) into Eq. (41), we obtain 


7r<n> = {-3 x±ufK 2 (x±) + X±G± (3x± ufK 3 (x±) ± 2e x± Cf) 


B 


± (8 (h±) 2 - 2) e~ x ±Cf - 15x1^/13 (x±) 


2e x± n±Ctuf n± n± , 

±- , ± 7T ± 


f3e~ x± Ctuf 

V X± K 3 (Xd 


X±I< 2 (x±) 2x± 2x± 

Substituting ft in Eq. (38) into Eq. (40), we obtain 


<n> 


(43) 


p-x±n ±~± 
e o 3 7T (11) 


1 


± 4..2 rt [ ( t j ± [^2 


f h±e x± (x± + 1) 


n±e 


-x± 


4x±^ 3 (Xd 

(x± + 1 ) 


2x1^2 (x±) 


2 (x±) -^2 (Xd 

9± 


± Jt {x ±e -fif (x± + 2 ) (Cf + 3)} + 


{xlCfK, (x±)) = 0. 

hi 


(44) 
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Finally, Fourier law, which is obtained from Eqs. (5), (27) and (28), yields the relation 


±\ 2 


0 ± 


I = —Kn 


(Bt) 


8 X ± (X± + 2) K 2 (2y±) + 5x^3 (2y±) ^/ a '0 ± + b r 


(45) 


From Eqs. (42)-(45), we obtain eight nonlinear equations for eight parameters (x±, , 7r^ 

and g±). Finally, we obtain following sets of eight parameters in cases (E) and (F) 


case (E) 

X- = 19.77, x+ = 2.15, q ] _ = -0.004371, q\ = -0.004222, 

7r ( “ n) = 0.007272, = -0.004014, uf = 0.007632, uf = 0.004955. 

case (F) 

X- = 97.15, x+ = 10.3, qi = -0.0004078, q\ = -0.0004243, 

7r ( - n) = 0.001399, 7fJ 1} = -0.0005426, uf = 0.002209, uf = 0.005246. 

From above data, we obtain (A6 1 ! 4 , Ad ^ 4 ) = (0.06712, 0.0186) in case (E), and ( Ad 14 , A^ 4 ) = 
(0.03886,0.02893) in case (F). As a result, (Ad} 4 , Ad^) in cases (E) and (F) are markedly 
different from that in tests (i) of cases (E) and (F). Provided that we can neglect no perme¬ 
ations in normal directions to two walls in Eq. (44), we obtain following sets of six parameters 
by setting (uf , ) = (—0.0121,0.0395) in case (E) and (uf,uf) = (—0.0082,0.046) in case 

(F) in Eqs. (42), (43) and (45) 


case (E) 

= 23.88, x+ = 2.407, q\ = -0.00369, q\ = -0.003632, 

7f ( - n) = -0.0005567, ft+ n) = -0.00448 
case (F) 

= 120.4, x+ = 13.28, q\ = -0.0002774, q\ = -0.0002873, 

7f ( " 1} = -0.00296, 7r+ n) = -0.001212. 

From above data, we obtain (Ad} 4 , AP) 4 ) = (—0.1156, 0.1232) in case (E), and (Ad 14 , A^ 4 ) = 
(—0.1614,0.247) in case (F), which are similar to (A6L, Ad + ) in tests (i) of cases (E) and 
(F), as shown in Tab. 1 . Meanwhile, we must remind that artificial sets of (uf,uf) violate 
no permeations in normal directions to two walls. Finally, we cannot describe temperature 
jumps on two walls in the framework of Grad’s 14 moments coupled with Fourier law in tests 
(i) of cases (E) and (F). Consequently, higher order moments beyond Grad’s 14 moments 
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must be considered to describe temperature jumps on two walls owing to the increase of f 
in tests (i) of cases (E) and (F). 

Substituting 6L and 9 +1 which are numerically obtained in tests (i) of cases (A)-(F), into 
Eqs. (25)-(28), p versus 1 — x\ and 9 versus 1 — x\ are plotted, as shown in Fig. 5. As 
shown in Fig. 5, profiles of p and 9 in tests (i) of cases (A)-(D) are similar to p and 9, which 
are calculated using Eqs. (25) and (26), namely, solutions of p and 9 under the thermally 
ultrarelativistic limit, whereas profiles of p and 9 in tests (i) of cases (E) and (F) are similar 
to p and 9, which are calculated using Eqs. (27) and (28), namely, solutions of p and 9 under 
the nonrclativistic limit. Consequently, we can conclude that the analytical solutions of the 
pure conductive state in Eqs. (25)-(28) are almost accurate, when we fix 9± to numerical 
values obtained using the relativistic Boltzmann equation. Therefore, the criticism for the 
use of the relativistic NSF equation in Eqs. (l)-(5) owing to its acausality and instability is 
presumably optimistic, when we discuss the conductive state of thermally relativistic fluids 
in flow-fields except for those in vicinities of two walls. Investigating profiles of 9 in detail, 
profiles of 9, which are numerically obtained in tests (i) of cases (A)-(C), are higher than the 
profile of 9 , which is calculated using Eqs. (25) and (26), whereas 9, which are numerically 
obtained in tests (i) of cases (A)-(C), increase more markedly toward the hot wall than 9, 
which is calculated using Eqs. (25) and (26). As a result, the approximation of q 1 on the 
basis of Fourier law in Eq. (24) is inadequate to demonstrate the profile of 9 around the hot 
wall owing to strongly nonequilibrium state around the hot wall. 

Temperature jumps of the nonrelativistic gas on two walls can be described using the first 
order theory on the basis of the NSF law, as discussed by Stefanov et al j^], when Kn = 0.01. 
Meanwhile, our numerical results of the conductive state indicate (A 9_, A 9 + ) = (0.24, 0.037) 
for Kn = 0.01 and (A6L, A6* + ) = (0.92,0.183) for Kn = 0.1 under the nonrelativistic limit, 
when 0+ = 1.0, 9~ = 0.1, Fr = oo and f ~ Kn/ \/~9. 

Here, we emphasize that temperature jumps on two walls under y -C 1 cannot be described 
in the framework of Grad’s 14 moments coupled with Fourier law even with modified trans¬ 


port coefficients proposed by Denied et al. llj, which are different from those obtained on 


the basis of Israel-Stewart theory. The transition of the conductive state to the convective 
state might be described in the framework of the relativistic NSF equation, when we set 
Kn as Kn <C 0.01, whereas the calculation using Kn <C 0.01 is markedly difficult because 
of the decrease of the time step in accordance with the decrease of Kn even with the most 
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advanced supercomputers. 


B. Numerical results of Rayleigh-Benard convection 

From above discussions on the conductive state, flow-fields in vicinities of two walls are 
under strongly nonequilibrium states, which do not allow us to describe the transition of the 
conductive state to the convective state of thermally relativistic fluids in the framework of 
the relativistic NSF equation. Actually, temperature jumps on two walls cannot be described 
in the framework of Grad’s 14 moments equation coupled with Fourier law, as discussed in 
Sec. IV-A. Therefore, the theoretical value of the critical Rayleigh number in appendix A, 
which is calculated by the relativistic NSF equation in Eqs. (l)-(5), cannot characterize the 
transition of the conductive state to the convective state with good accuracy. We, however, 
use the definition of the Rayleigh number (Ra) to characterize such a transition, because 
we cannot introduce a new parameter other than Ra, which characterizes such a transition 
with better accuracy than Ra. 

The critical Rayleigh number (Ra cr ), which corresponds to Ra in tests (iii), decreases, as x~ 
(x+) decreases, when 0.1 < x~ A 10, whereas Ra cr decreases, as x~ (x+) decreases, when 
2.9 x 10~ 3 < X- < 0.1, as shown in Tab. 1. The critical Fr (Fr cr ), which corresponds to Fr 
in tests (iii), is in the range of 30.36 < Fr cr < 36.22, whereas Fr cl . ~ 30 under nonrelativistic 
limit, as obtained by Stefanov et al Q. 

A6L increases, as g increases in cases (A), (B), (D), (E) and (F), whereas A6L under g = 5.4 
(test (iii)) is larger than A6L under g = 5.5 (test (iv)), as shown in Tab. 1. A 9 + under 
the conductive state (test (i)) is larger than A 6 + under the weakly wavy state (test (ii)) in 
cases (A)-(E), whereas A 6 + under the conductive state (test (i)) is slightly smaller than A 6+ 
under the weakly wavy state (test (ii)) in case (F), whereas A0+ 4 under the pure conductive 
state increases, as g increases, as shown in Tab. 1. Finally, A 6 + increases, as g increases 
from the weakly wavy state. 


V. CONCLUDING REMARKS 

We investigated the thermal conduction and convection of thermally relativistic fluids 
under the gravitational force between two parallel walls, theoretically and numerically, when 
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the Lorentz contraction of thermally relativistic fluids is negligible. Numerical results of 
the thermal conduction of thermally relativistic fluids under the gravitational force cannot 
be described in the framework of Grad’s 14 moments coupled with Fourier law owing to 
strongly nonequilibrium states on the hot and cold walls, when the flow-held is under the 
transition regime between the rarefied and continuum regimes. Additionally, absolute values 
of temperature jumps on two walls increase, as the temperature decreases in the range 
of 2.9 x ICG 4 < x < 120, because the representative relaxational time increases as the 
temperature decreases. Meanwhile, profiles of the density and temperature between two 
walls under the conductive state shows good agreements with their analytical solutions, 
which are obtained using the relativistic NSF equation, when we fix temperatures on two 
walls to those, which are numerically obtained. Therefore, we conclude that how-fields except 
for those in vicinities of two walls can be described using the relativistic NSF equation owing 
to their weakly nonequilibrium states with good accuracy. As a result of such strongly 
nonequilibrium states in vicinities of two walls beyond the relativistic NSF equation, the 
Rayleigh number cannot be a parameter, which characterizes the transition of the conductive 
state to the convective state with good accuracy, when the how-held is under the transition 
regime between the rarehed and continuum regimes, whereas the Rayleigh number can be 
a parameter, which characterizes the transition of the conductive state to the convective 
state with good accuracy under the nonrelativistic limit, when the how-held is under the 
transition regime between the rarehed and continuum regimes. 


Appendix A: Critical Rayleigh number for thermally relativistic fluids 


The temperature prohlc is expressed with 


9 (xi) = 9 0 + xi (0i - 6>o) + 59 (xi,x 2 ) (9 1 - 9 0 ) (0 < x 1 < 1), (Al) 

'-«v-' 

<50 


where 9q = 9- and 9\ — 9 + . 

To eliminate the gradient of the pressure in Eq. (20), we rewrite Eq. (20) using the vorticity. 
Equation (20) can be rewritten as 


1 

Pr 


85 u 
~df 


5u ■ 5u 


5<jj x 5u 


= -V5p' - W 2 Su + Ra 59e, 


(A2) 
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where we used the relation bu ■ \7bu = V ( Su 9 Su ) + bu> x bu. 

Taking the curl of Eq. (A2) 

1 /<9V x bu \ 

^— + Vxtox^u = V 2 (V x 5u) + Ra (V x S9e) , (A3) 

Pr \ dt J 

where we used the relation that the curl of a gradient is equal to zero. 

Substituting buj = V x bu and V x bu: x bu = bu ■ Vbuj — buj ■ X7bu into Eq. (A3), we 
obtain 


1 / dbtv 


Pr V dt 

Two dimensional flow yields following relations 

. d§U 2 dbu\ 

buj 3 = 


bu ■ Vbu> - bu: ■ Vbuj = V 2 to + Ra (V x bde) , 


(A4) 


dx\ dx2 

bCj ■ Vbu = 0, 

dbd 


e 3 , bui = bu 2 = 0, 


V x bde = 


dx< 


e 3 , (ve = (1,0,0)) 


Equation (A4) is rewritten using Eq. (A5) as 

1 / dbti 3 


+ bu ■ V(fu; 3 ) = V 2 bu 3 — Ra 


dbd 


Pr V 9t J dx 2 

Substituting Eq. (14) and = XL into Eq. (13), we obtain 

db6 r _ c pr r rt 

bu ■ X7b9 = \7~b9, 


dt 


(AS) 


(A6) 


(AT) 


Substituting Eq. (Al) into Eqs. (A6) and (A7), we obtain by neglecting all the nonlinear 
terms 


1 dbu 3 ^ dbd' 

T — = V &3 - Ra^-, 


dt 


bui = VW, 


bui, bu 2 and bco 3 are expressed with the stream function 0 as 

d ^ r- ^0 r n2 / 

bu\ = , bu 2 = bui 3 = -V 0. 

OX 2 OX\ 

Substituting Eq. (A10) into Eqs. (A8) and (A9), we obtain 

1 9V 2 0 ,^2 a T. dSff 

= V 2 (V 2 0) + Ra 


p r dt 
dbd' 90 
dt dx 2 


dx 2 ’ 


= VW. 


(AS) 

(A9) 

(A10) 

(All) 

(A12) 




















We set 0 and 50' as 
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0 = 0 (xi) exp (at — 1 KX 2 ), 
50' = 50' (xi) exp (at — ikx 2 ). 


Substituting Eqs. (A13) and (A14) into Eqs. (All) and (A12), we obtain 

2 


p r a \dx\ 


a ( —n — k 2 ) 0 = —Rai/c50' + ( -Ay — k 2 ) 0, 


\dxf 


<t 50' — i/c0 = — ( ^ / - 


cr > 0 yields the instability and a < 0 yields the stability of the convection, 
consider the case of a = 0. 

Substituting a = 0 into Eqs. (A15) and (A16), we obtain 

K ^ (s,)= s;(4r' c2 ) 

Substituting 0 (ah) = exp (£xi) into Eq. (A17), we obtain six roots 


fo,i = ±*« j 

f* 

6,3 = | 

A 

£4,5 = ±ire j 

[■♦< 


VRa J 

( AC 2 \ ® 


K 




exp 


7T 


From Eq. (A18), 0 = 0 A* exp (&X0 (i = 0,1, 2, 3,4, 5). 

We consider three types of the boundary condition as follows 


Rigid boundary for two walls 

uf = 0, 0(0) = 0_, 0(1) = 0+. 

Slip boundary for two walls 
dTi^~ 

_A = 0, 0(0) = 0_, 0(1) = 0 + . 

dx\ 

Slip boundary for hot wall 

/7t7+ 

P- = 0. aj = 0, 0(0) = o_, o(i) = o+. 


(A13) 

(A14) 

(A15) 
(A16) 
Then, we 

(A17) 


(A18) 

(A19) 

(A20) 


(A21) 
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The case of the slip boundary for the hot wall and rigid boundary for the cold wall might 
be significant, when the density on the hot wall is rarefied by the transfer of fluids from the 
hot wall to cold wall. 

From Eqs. (A15), (A16), (A19), (A20) and (A21), we can rewrite Eqs. (A19)-(A21) as 

Equation (A19) —> — 0, = 0, ( — k 2 ^) i>± = 0, (A22) 

dx i y dx ^ J 

(fib, „ { d 2 \ 2 - 

Equation (A20) —> 9 = 0, i/j± = 0, ( —j — k 2 ) z/q_ = 0, (A23) 

dx y dx J 

Equation (A21) —> ^ 9 + = 0, -j— = 0, ijj± = 0, (—r — k 2 ^ = 0. (A24) 

dx ^ dx x y dx ^ J 

Substituting ^ = X^=o A ex P (£i^i) (® — 0,1, 2, 3,4, 5) into Eqs. (A22)-(A24), we obtain 

M A = 0, (A25) 

where A = t (A i ) (i = 0,1, 2, 3,4, 5). 

Finally, det M = 0 yields the relation between Ra and k, as shown in Fig. 6. Finally, we 
obtain Ra cr for three types of the boundary condition from Fig. 6, as 


Rigid boundary for two walls 

Ra cr = 1708. (A26) 

Slip boundary for two walls 

Ra cr = 657. (A27) 

Slip boundary for hot wall and rigid boundary for cold wall 

Ra cr = 1100. (A28) 

These values of Ra cr for three types of the boundary condition are quite same as those of 
nonrelativistic fluids. 


Appendix B: Numerical method of solving relativistic Boltzmann equation 

The Relativistic Boltzmann equation (RBE) is written as jlo| 

df 


P a i6;=Q(fJ) 


dx 


[ [ (/:/' - /./) iv<M— 

Im 3 Jn 2 P *o 


(Bl) 
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where / = / (t,x l ,p l ) (i = 1,2,3), and F is the Lorentz invariant flux. In Eq. (Bl), terms 
with a prime indicate conditions after collisions, f 2 2 is the solid angle space and & 3 is the 
momentum space stretched by {^ 3 |(—oo, — oo, — oo) < (p 1 ,^ 2 ,^ 3 ) < (oo, oo, oo)}. 
x a , p a and F are given by 


x a =[ct, x 1 , x 2 , x 3 ) , 
p a —m n f (v) (c, v 1 , v 2 ,u 3 ), 


(B2) 

(B3) 


p°pl p°pl 

F= - g 0 = - 


v - u *) 2 - — (v x v *) 2 


g 0 =\l ( v - u *) 2 - — (v x v ^ 2 


=\/ ( P*Pa f - rri 4 c 4 . 


(B4) 


In Eq. (B3), 7 (f) is the Lorentz factor, which is defined by 7 (v) = 1/y/l — v 2 /c 2 . c is 
the speed of light and v l (i = 1,2,3) is the ith component of the particle velocity vector 
v = (n 1 ,n 2 ,n 3 ). I 11 Eqs. (B3) and (B4), m is the molecular mass. In Eq. (B4), g 0 is Mpller’s 
relative velocity. In Eq. (Bl), a is the differential cross section and dfl is the solid angle 
element. In Eqs. (Bl) and (B4), terms with an asterisk subscript belong to the collision 
partner. 

Rewriting Eq. (Bl) in Lorentz variant form yields 


df , i df 

-h V -r 

dt dx l 


/^3 Jn 2 


(/;/'- /./) gt<rdSld 3 p, 


(B5) 


The direct simulation Monte Carlo (DSMC) method 9] developed for nonrelativistic gases 
can be extended into relativistic gases. The Courant-Friedlichs-Lewy (CFL) condition 12], 
which is required on the left-hand side of Eq. (Bl), requires that the time step At approxi¬ 
mates to zero when p a —> 00 , namely v —> c. Consequently, the left-hand side of Eq. (Bl) 
leads to the numerical stiffness via At —y 0 when we consider thermally relativistic flow. 
Thus, Eq. (B5) is solved instead of Eq. (Bl). As a numeric of the collision term in Eq. 
(B5), the majorant frequency scheme 13J is used, whereas the past work 14] 15] used the 
Bird’s scheme jg], which calculates the collision-probability for all binary collision-pairs in 
the numerical cell. Consequently, the computational time required by the Bird’s scheme is 
markedly longer than that by required by the majorant frequency scheme, because the ma¬ 
jorant frequency scheme calculates the collision-pair the maximum collision number times. 
In the majorant frequency scheme, the maximum collision number during At is obtained for 
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a hard-sphere particle from Eq. (5) as 

*W = ^ (N - 1) np(u)a T (c/ 0 ) max At = (N - 1) ncp(u)a T At (v (c/ 0 ) max = 2c), (B6) 

where iV is the number of sample particles in the cell. A collision pair is selected i/ max 
times. The two particles selected induce a binary collision when the random number W 
(0 < W < 1) satisfies 

—^- = — < yT, (B7) 

(# 0 ) max 2 C 

where g 0 is the Mqller’s relative velocity for two particles selected. 

Before and after the binary collision between particles 1 and 2, the total energy and total 
momentum of the binary particles must be conserved. The conservation of the energy E 
and the momentum p = ( p 1 ,p 2 ,p 3 ) before and after a binary collision is written as 


E + E*—E' + E'* — E tot , 

(B8) 

P + P*=P + P* = Ptot, 

(B9) 


where E = my(r0c 2 , E * = my (u*) c 2 , p = my (T) v and p* = my (u*) u*. 

E and p are related as follows 

E = . (BIO) 

In this paper, a binary collision is calculated using the following algorithm: 

(a) Calculate the total energy (E tot ) and total momentum (p tot = (p' <)t , p 2 )t , p 3 ot )) of binary 
particles from the left-hand sides of eqs. (B8) and (B9). 

(b) Redistribute the energy to binary particles, namely, E' and A', using a random number 
on the right-hand side of Eq. (B8). From eq. (BIO), the norms of momenta, namely, \p'\ 
and |p'J, are fixed. 

(c) Decide the direction of the momentum by the law of cosines so that the total momentum 
is conserved in Eq. (B9) as 


^ P l ' ^ 


0 1 


( P* ^ 


0 1 

P 2 ' 

= M(<j>,tp)N(r) 

— \p'\ sin 0i 


P* 

— M (0, ip) N (r) 

— \P* \ sin 02 

V ) 


^ \p'\ COS 0i y 




V |p* | COS 02 y 


(Bll) 
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where 0i, 02 , M (0, <p) and N ( r ) are given by 


= arccos 


Ip' 

l 2 + 

Ptot 

2 

lp*l 

n 


2 1 

p'l 

IPtot 



= arccos 


I / |2 , I 1 2 

If* I + iPtot 


Ip 


2 |p* I |p t , 


(B12) 



COS 0 COS (p 

— sin tp 

sin 0 cos ip 


cos r 

— sin r 

0 

M (0, ip) — 

cos 0 sin p 

cos p 

sin 0 sin ip 

, N(r) = 

sin r 

cos r 

0 


, — sin 0 

0 

COS 0 , 


0 

0 

1 


, (B13) 


where r is a random number in the range of 0 < r < 27T, and 0 and ip are given by 


0 = arccos 


pL 

IPtot 


(B14) 


if 0 < = arccos 


Ptot 


Ipt, 


sm < 


else 27 t — arccos ( T —^. (B15) 

\ | Plot | sm 0 J 

In our numerical analysis, 1638400 sample particles are simulated to solve the RBE in the 
numerical domain in Fig. 1 using 128 processors (4 nodes) of HITACHI SR 16000 through 
the collaborative work with Information Technology center of University of Tokyo. 


Appendix C: Derivation of relativistic Navier-Stokes-Fourier equation in Eqs. (l)-(3) 


Balance equations of the mass, momentum density, and energy density are written in 


Eckart’s frame as 


iq| 


Dn + nV a U a = 0, 



DU a 


V“ (p + n) - Vp7 r< Q/J) 


- n DU a 


—Dq a - q a VpUP - q^VpU a - ^U a q^DU p - U a n { ^ ] VpUA , 

nDe = - (p + n) + 7 T {ap) VpU a - V a q a + 4 q a DU a , 

& 


(Cl) 


(C2) 

(C3) 


where n = p/m is the number density, U a = y(u) (c, v/) (i— 1,2,3, 7 (v) = 1 /yjl — v 2 /c 2 : 
Lorentz factor) is the four flow velocity, D = U"V cv is the convective time derivative, 
V" = A 0 / 3 03 , in which = (r— U a U^ /c 2 ) 9^, where rj a P = diag (1,-1,—1,-1), is the 
projector, e is the energy density, and h E = mc 2 G is the enthalpy per a particle. 

In this paper, we investigate thermal fluctuations under static state in the laboratory frame. 
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Then, we assume that the product of nonlinear terms, which are expressed by products of U a 
(or U a ) and II or q a in Eqs. (C2) and (C3), are negligible owing to Vi <C c (i = 1, 2, 3). 
Consequently, linearized balance equations of the momentum density and energy density are 
written from Eqs. (C2) and (C3) by neglecting nonlinear terms as 10] 


nhi 


■DU a = V“ (p + n) - V pTT 


(ap) 


jDq a , 


nDe = -pV p U p - VpqP. 


(C4) 

(C5) 


In Eq. (C4), we assume that th term — Dq a /c 2 is negligible, and rewrite Eq. (C4) as 

^ DU a = V Q (p + n) - Vfl. (C6) 

er 

From Eqs. (Cl), (C5) and (C6), we readily obtain Eqs. (l)-(3) using relations, De = c v D9, 
Jie = mc 2 G, and 7 (v) = 1. 
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Figure captions 

Figure 1: Schematic of flow-field. 

Figure 2: A6 M versus x~ and A 9^ versus x+ obtained using Eq. (36), when g = 10 1 , 100 
and 1000 (left frame). AO 14 versus x~ and A 9^ versus x+ obtained using Eq. (37), when 
g = 10 -4 , 10^ 3 and 10 -2 (right frame). 

Figure 3: U\ versus 1 — X\ in y v axis and II in y 2 axis versus 1 — X\ in tests (i) of cases 
(A)-(F). 

Figure 4: q 1 , g 4 SP , fr^n), and versus 1 — X\ in tests (i) of cases (A)-(F). 

Figure 5: 9 versus 1 — X\ in y\ axis and p versus l — x\ in y 2 axis, where analytical solutions 
of 6 and p in Eqs. (25)-(28) are calculated using and 9 + , which are numerically obtained. 
Figure 6: Ra versus n for three types of the boundary condition, where shaded areas corre¬ 
spond to instable regimes. 
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TABLE I. Numerical results of cases (A)-(F) 


Case (A) 

9 


Ad + 

A0 U (Eq. (36)) A9 

+ (Eq- (36)) 

Ra 

Fr 

Status 

test (i) 

100 

0.059 

0.01442 

0.03352 

0.02605 

437.7 

68 

C 

test (ii) 

193.5 

0.06285 

0.0137 

0.03352 

0.02605 

846.1 

35.15 

ww 

test (iii) 

194 

0.0681 

0.01485 

0.03352 

0.02605 

848.6 

35.02 

4-V 1 

test (iv) 

195 

0.075 

0.02 

0.03352 

0.02605 

856.2 

34.66 

4-V 2 

Case (B) 

9 

A 

a e + 

Adi 4 (Eq. (36)) Ad 

+ (Eq- (36)) 

Ra 

Fr 

Status 

test (i) 

0 

0.064 

0.0155 

0.05768 

0.02579 

0 

oo 

C 

test (ii) 

19.15 

0.068 

0.01222 

0.05766 

0.0258 

835.7 

35.57 

WW 

test (iii) 

19.25 

0.07 

0.01271 

0.05766 

0.0258 

840.3 

35.38 

4-V 1 

test (iv) 

19.5 

0.076 

0.017 

0.05766 

0.0258 

853.8 

34.76 

4-V 2 

Case (C) 

9 

A0„ 

Ad + 

Adi 4 (Eq. (36)) A9 

+ (Eq- (36)) 

Ra 

Fr 

Status 

test (i) 

1 

0.06 

0.022 

0.1264 

0.02508 

151.9 

195.6 

C 

test (ii) 

5.35 

0.0665 

0.02031 

0.1264 

0.0251 

811.1 

36.62 

WW 

test (iii) 

5.4 

0.07143 

0.02186 

0.1264 

0.0251 

819.3 

36.22 

4-V 1 

test (iv) 

5.5 

0.069 

0.0222 

0.1264 

0.0251 

835 

35.6 

4-V 2 

Case (D) 

9 

A 

a e + 

Adi 4 (Eq. (37)) Ad 

+ (Eq- (37)) 

Ra 

Fr 

Status 

test (i) 

0.1 

-0.01927 

0.0336 

-0.02285 

0.02309 

381 

77.31 

C 

test (ii) 

0.235 

-0.01674 

0.033 

-0.02244 

0.0235 

894.5 

32.92 

WW 

test (iii) 

0.24 

-0.01252 

0.0337 

-0.02243 

0.02351 

913.6 

32.21 

4-V 1 

test (iv) 

0.25 

0.00014 

0.036614 

-0.0224 

0.02355 

952.4 

30.83 

4-V 2 

Case (E) 

9 

A 0_ 

a e + 

Adi 4 (Eq. (37)) A9 

+ (Eq- (37)) 

Ra 

Fr 

Status 

test (i) 

0.01 

-0.12 

0.121 

-0.000234 

-0.00954 

395.1 

83.3 

C 

test (ii) 

0.025 

-0.11 

0.105 

-0.0002054 

-0.01106 

969.6 

33.92 

WW 

test (iii) 

0.0275 

-0.11 

0.1193 

-0.0002012 

-0.01135 

1083 

30.36 

4-V 1 

test (iv) 

0.03 

-0.076 

0.1248 

-0.000197 

-0.01164 

1184 

27.65 

4-V 2 

Case (F) 

9 

A 

a e + 

Adi 4 (Eq. (37)) A9 

+ (Eq- (37)) 

Ra 

Fr 

Status 

test (i) 

0.0025 

-0.16 

0.2494 

-6.878 xlO" 6 

-0.0006361 

608.24 60.05 

C 

test (ii) 

0.0045 

-0.159 

0.25 

-4.984xl0 -6 

-0.001073 

1095 

33.33 

WW 

test (iii) 

0.00475 

-0.15 

0.2513 

-4.805 xlO" 6 

-0.001149 

1290 

31.52 

4-V 1 

test (iv) 

0.005 

-0.14 

0.2527 

-4.636 xlO" 6 

-0.001232 

1218 

29.9 

4-V 2 
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